/*
 * SelfComputeIsotropicElasticityTensor.C
 *
 *  Created on: 2018年1月9日
 *      Author: liuxiao
 */
/****************************************************************/
/* MOOSE - Multiphysics Object Oriented Simulation Environment  */
/*                                                              */
/*          All contents are licensed under LGPL V2.1           */
/*             See LICENSE for full restrictions                */
/****************************************************************/
#include "SelfComputeIsotropicElasticityTensor.h"
registerMooseObject("TrilobitaApp", SelfComputeIsotropicElasticityTensor);
template<>
InputParameters validParams<SelfComputeIsotropicElasticityTensor>()
{
  InputParameters params = validParams<ComputeElasticityTensorBase>();
  params.addClassDescription("Compute an isotropic elasticity tensor.");
  params.addRequiredCoupledVar("temperature", "Coupled Temperature");
  params.addParam<std::vector<Real>>("bulk_modulus_list", "The bulk modulus vector for the material.");
  params.addParam<std::vector<Real>>("lambda_list", "Lame's first constant vector for the material.");
  params.addParam<std::vector<Real>>("poissons_ratio_list", "Poisson's ratio vector for the material.");
  params.addParam<std::vector<Real>>("shear_modulus_list", "The shear modulus vector of the material.");
  params.addParam<std::vector<Real>>("youngs_modulus_list", "Young's modulus vector of the material.");
  params.addParam<std::vector<Real>>("bulk_modulus_T_list", "The temperature list to compute bulk modulus for the material.");
  params.addParam<std::vector<Real>>("lambda_T_list", "The temperature list to compute Lame's first constant for the material.");
  params.addParam<std::vector<Real>>("poissons_ratio_T_list", "The temperature list to compute Poisson's ratio for the material.");
  params.addParam<std::vector<Real>>("shear_modulus_T_list", "The temperature list to compute The shear modulus of the material.");
  params.addParam<std::vector<Real>>("youngs_modulus_T_list", "The temperature list to compute Young's modulus of the material.");
  return params;
}

SelfComputeIsotropicElasticityTensor::SelfComputeIsotropicElasticityTensor(const InputParameters & parameters) :
    ComputeElasticityTensorBase(parameters),

    _bulk_modulus_set( parameters.isParamValid("bulk_modulus_list") ),
    _lambda_set( parameters.isParamValid("lambda_list") ),
    _poissons_ratio_set( parameters.isParamValid("poissons_ratio_list") ),
    _shear_modulus_set( parameters.isParamValid("shear_modulus_list") ),
    _youngs_modulus_set( parameters.isParamValid("youngs_modulus_list") ),
	 _temperature(coupledValue("temperature")),
	_d((0.0,0.0)),
	_bulk_modulus_list( _bulk_modulus_set ? getParam<std::vector<Real>>("bulk_modulus_list"):_d),
	_bulk_modulus_T_list( _bulk_modulus_set ? getParam<std::vector<Real>>("bulk_modulus_T_list"):_d),
	_lambda_list( _lambda_set ? getParam<std::vector<Real>>("lambda_list") :_d),
	_lambda_T_list( _lambda_set ? getParam<std::vector<Real>>("lambda_T_list") : _d),
	_poissons_ratio_list( _poissons_ratio_set ?  getParam<std::vector<Real>>("poissons_ratio_list") :_d),
	_poissons_ratio_T_list( _poissons_ratio_set ?  getParam<std::vector<Real>>("poissons_ratio_T_list") : _d),
	_shear_modulus_list( _shear_modulus_set ? getParam<std::vector<Real>>("shear_modulus_list") :_d),
	_shear_modulus_T_list( _shear_modulus_set ? getParam<std::vector<Real>>("shear_modulus_T_list") :_d),
	_youngs_modulus_list( _youngs_modulus_set ? getParam<std::vector<Real>>("youngs_modulus_list") : _d),
	_youngs_modulus_T_list( _youngs_modulus_set ? getParam<std::vector<Real>>("youngs_modulus_T_list") :_d),
	_bulk_modulus(declareProperty<Real>("bulk_modulus")),
    _lambda(declareProperty<Real>("bulk_modulus")),
	_poissons_ratio(declareProperty<Real>("poissons_ratio")),
	_shear_modulus(declareProperty<Real>("shear_modulus")),
	_youngs_modulus(declareProperty<Real>("youngs_modulus"))
{

	if (_shear_modulus_set && _bulk_modulus_set)
	{
		_func_bulk_modulus_T =  LinearInterpolation(_bulk_modulus_T_list, _bulk_modulus_list);
//		_func_lambda_T =  LinearInterpolation(_shear_modulus_T_list, _shear_modulus_list);
//		_func_poissons_ratio_T =  LinearInterpolation(_shear_modulus_T_list, _shear_modulus_list);
		_func_shear_modulus_T =  LinearInterpolation(_shear_modulus_T_list, _shear_modulus_list);
//		_func_youngs_modulus_T =  LinearInterpolation(_youngs_modulus_T_list, _youngs_modulus_list);
	}
	else if (_lambda_set && _shear_modulus_set)
	{
//		_func_bulk_modulus_T =  LinearInterpolation(_bulk_modulus_T_list, _bulk_modulus_list);
		_func_lambda_T =  LinearInterpolation(_lambda_T_list, _lambda_list);
//		_func_poissons_ratio_T =  LinearInterpolation(_shear_modulus_T_list, _shear_modulus_list);
		_func_shear_modulus_T =  LinearInterpolation(_shear_modulus_T_list, _shear_modulus_list);
//		_func_youngs_modulus_T =  LinearInterpolation(_youngs_modulus_T_list, _youngs_modulus_list);
	}
	else if(_youngs_modulus_set && _poissons_ratio_set)
	{
		_func_poissons_ratio_T =  LinearInterpolation(_poissons_ratio_T_list, _poissons_ratio_list);
		_func_youngs_modulus_T =  LinearInterpolation(_youngs_modulus_T_list, _youngs_modulus_list);
	}
}

void
SelfComputeIsotropicElasticityTensor::computeQpElasticityTensor()
{

	  std::vector<Real> iso_const(2);
	  if (_lambda_set && _shear_modulus_set)
	  {
		  _lambda[_qp]=_func_lambda_T.sample(_temperature[_qp]);
		  _shear_modulus[_qp]=_func_shear_modulus_T.sample(_temperature[_qp]);
		  iso_const[0] = _lambda[_qp];
		  iso_const[1] = _shear_modulus[_qp];
	  }
	  else if (_youngs_modulus_set && _poissons_ratio_set)
	  {

		  _youngs_modulus[_qp]=_func_youngs_modulus_T.sample(_temperature[_qp]);
		  _poissons_ratio[_qp]=_func_poissons_ratio_T.sample(_temperature[_qp]);
		  iso_const[0] = _youngs_modulus[_qp]*_poissons_ratio[_qp]/((1 + _poissons_ratio[_qp])*(1 - 2*_poissons_ratio[_qp]));
		  iso_const[1] = _youngs_modulus[_qp]/(2*(1 + _poissons_ratio[_qp]));
	  }
	  else if (_shear_modulus_set && _bulk_modulus_set)
	  {
		  _shear_modulus[_qp]=_func_shear_modulus_T.sample(_temperature[_qp]);
		  _bulk_modulus[_qp]=_func_bulk_modulus_T.sample(_temperature[_qp]);
		  iso_const[0] = _bulk_modulus[_qp] - 2.0/3.0*_shear_modulus[_qp];
		  iso_const[1] = _shear_modulus[_qp];
	  }
	  else
	    mooseError("Incorrect combination of elastic properties in ComputeIsotropicElasticityTensor. Possible combinations are: lambda and shear_modulus, youngs_modulus and poissons_ratio, or bulk_modulus and shear_modulus.");

	  //Fill elasticity tensor
	  _Cijkl.fillFromInputVector(iso_const, RankFourTensor::symmetric_isotropic);
  //Assign elasticity tensor at a given quad point
  _elasticity_tensor[_qp] = _Cijkl;
}





